
##########################################################################################################################
#### Replication Data for Figures in "Ethnicity and Response to Internal Environmental Migrants in the United States"#####
#### Forthcoming in Political Research Quarterly #########################################################################
##########################################################################################################################


#setwd("C:/Users/kash8/OneDrive/Documents/Climate Stress Response Project")
setwd("~/ash_winchester_prq_replication")

library(ggplot2)
library(stringr)

#################################
###ARTICLE FIGURES###############
#################################

####Figure 2#####


q30<- read.csv("q30.csv",sep = "\t")


White <- data.frame(Treatment = q30$race,
                          Difference = q30$treat_white,
                          SE = q30$se_white,
                          modelName = "White") 

Black <- data.frame(Treatment = q30$race,
                          Difference = q30$treat_black,
                          SE = q30$se_black,
                          modelName = "Black") 
Latinx <- data.frame(Treatment = q30$race,
                          Difference = q30$treat_latinx,
                          SE = q30$se_latinx,
                          modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Approval of Migrants from Control") + labs(color='Treatment')
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure 3#####


ethnocentrism<- read.csv("ethnocentrism.csv", sep="\t")


White <- data.frame(Treatment = ethnocentrism$race,
                    Difference = ethnocentrism$treat_white,
                    SE = ethnocentrism$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = ethnocentrism$race,
                    Difference = ethnocentrism$treat_black,
                    SE = ethnocentrism$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = ethnocentrism$race,
                     Difference = ethnocentrism$treat_latinx,
                     SE = ethnocentrism$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Ethnocentrism from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure 4####


tolerance<- read.csv("tolerance.csv", sep="\t")


White <- data.frame(Treatment = tolerance$race,
                    Difference = tolerance$treat_white,
                    SE = tolerance$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = tolerance$race,
                    Difference = tolerance$treat_black,
                    SE = tolerance$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = tolerance$race,
                     Difference = tolerance$treat_latinx,
                     SE = tolerance$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Tolerance from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####################################
####Figure 5: Intergroup Contact####
####################################




###Changes in perceived income and Education based on intergroup contact line####

conf=.95
title="Black locals' perceived income of white migrants"
xlabel="Percentage Migrants' Ethnicity within 5km of locals"
ylabel="Change in Perceived Income"


# Get coefficients of variables
beta_1 =  .9852312
beta_3 =   -2.260816   

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1.05



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =   .09733762  + (x_2^2)* .36494442     + 2*x_2* -.15491349 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("0%","100%"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)




###Changes in perceived income and Education based on intergroup contact line####

conf=.95
title="White locals' perceived income of Black migrants"
xlabel="Percentage Migrants' Ethnicity within 5km of locals"
ylabel="Change in Perceived Income"


# Get coefficients of variables
beta_1 =  -.3662998
beta_3 =  .8993736    

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1.05



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =       .01401458  + (x_2^2)* .17376841      + 2*x_2*-.02277798 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("0%","100%"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



conf=.95
title="White locals' perceived education of Black migrants"
xlabel="Percentage Migrants' Ethnicity within 5km of locals"
ylabel="Change in Perceived Education"


# Get coefficients of variables
beta_1 =  -.4392775
beta_3 =  .9339317     

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1.05



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .01430189  + (x_2^2)* .17639415     + 2*x_2*-.02309068 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("0%","100%"))

# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



conf=.95
title="Black locals' perceived income of Latinx migrants"
xlabel="Percentage Migrants' Ethnicity within 5km of locals"
ylabel="Change in Perceived Income"


# Get coefficients of variables
beta_1 =   1.061894
beta_3 =    -1.69461   

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1.05



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .05037156  + (x_2^2)*    .71155833     + 2*x_2*-.11665457 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("0%","100%"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



conf=.95
title="Black locals' perceived education of Latinx migrants"
xlabel="Percentage Migrants' Ethnicity within 5km of locals"
ylabel="Change in Perceived Education"

# Get coefficients of variables
beta_1 =   .8738795
beta_3 =   -2.175956    

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1.05



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .0497956  + (x_2^2)* .7497171     + 2*x_2*-.11808155 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("0%","100%"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)


#################################
###APPENDIX FIGURES##############
#################################

####Figure A1#####


q30<- read.csv("q30_max.csv",sep="\t")


White <- data.frame(Treatment = q30$race,
                    Difference = q30$treat_white,
                    SE = q30$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = q30$race,
                    Difference = q30$treat_black,
                    SE = q30$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = q30$race,
                     Difference = q30$treat_latinx,
                     SE = q30$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Approval of Migrants from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure A2#####

q32<- read.csv("q32_max.csv", sep="\t")


White <- data.frame(Treatment = q32$race,
                    Difference = q32$treat_white,
                    SE = q32$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = q32$race,
                    Difference = q32$treat_black,
                    SE = q32$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = q32$race,
                     Difference = q32$treat_latinx,
                     SE = q32$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Volunteer Time from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure A3#####

ethnocentrism<- read.csv("ethnocentrism_max.csv", sep="\t")


White <- data.frame(Treatment = ethnocentrism$race,
                    Difference = ethnocentrism$treat_white,
                    SE = ethnocentrism$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = ethnocentrism$race,
                    Difference = ethnocentrism$treat_black,
                    SE = ethnocentrism$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = ethnocentrism$race,
                     Difference = ethnocentrism$treat_latinx,
                     SE = ethnocentrism$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Ethnocentrism from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure A4#####

tolerance<- read.csv("tolerance_max.csv", sep="\t")


White <- data.frame(Treatment = tolerance$race,
                    Difference = tolerance$treat_white,
                    SE = tolerance$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = tolerance$race,
                    Difference = tolerance$treat_black,
                    SE = tolerance$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = tolerance$race,
                     Difference = tolerance$treat_latinx,
                     SE = tolerance$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Tolerance from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


####Figure A7#####

q32<- read.csv("q32.csv",sep = "\t")


White <- data.frame(Treatment = q32$race,
                    Difference = q32$treat_white,
                    SE = q32$se_white,
                    modelName = "White") 

Black <- data.frame(Treatment = q32$race,
                    Difference = q32$treat_black,
                    SE = q32$se_black,
                    modelName = "Black") 
Latinx <- data.frame(Treatment = q32$race,
                     Difference = q32$treat_latinx,
                     SE = q32$se_latinx,
                     modelName = "Latinx") 
FullModelFrame <- data.frame(rbind(White, Black, Latinx))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Treatment, ymin = Difference - SE*interval1,
                                ymax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Treatment, y = Difference, ymin = Difference - SE*interval2,
                                 ymax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("Respondent Race") + ylab("Change in Volunteer Time from Control") + labs(color='Treatment')
scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().






#####Figure A8: FIGURES FOR MEXICAN NATIONAL HERITAGE MODERATION####

conf=.95
title="Volunteer Time for White Migrants by Latinx Respondents"
xlabel=""
ylabel="Change in Willingness to Volunteer Time"
factor_labels=c("No Mexican Heritage","Mexican Heritage")


# Get coefficients of variables
beta_1 =   -.5832068  
beta_3 =    .6169158

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .03967843   + (x_2^2)*   .06824965     + 2*x_2*   -.01876497 

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)

conf=.95
title="White Migrants' Perceived Education by Latinx Respondents"
xlabel=""
ylabel="Change in Perceived Education"
factor_labels=c("No Mexican Heritage","Mexican Heritage")


# Get coefficients of variables
beta_1 =   -.1219883  
beta_3 =    .6003621

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .04030238  + (x_2^2)*  .06978449    + 2*x_2*   -.04038492

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



conf=.95
title="Black Migrants' Perceived Education by Latinx Respondents"
xlabel=""
ylabel="Change in Perceived Education"
factor_labels=c("No Mexican Heritage","Mexican Heritage")


# Get coefficients of variables
beta_1 =   -.0067674   
beta_3 =   .5019253 

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .03755447   + (x_2^2)*  .06748209   + 2*x_2* -.03755169

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)


conf=.95
title="Black Migrants' Perceived Employment by Latinx Respondents"
xlabel=""
ylabel="Change in Perceived Employment"
factor_labels=c("No Mexican Heritage","Mexican Heritage")


# Get coefficients of variables
beta_1 =    .245279    
beta_3 =   .5620891 

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =      .03773706   + (x_2^2)*  .06965061    + 2*x_2* -.03755152

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



